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Abstract 



We propose a high resolution finite volume scheme for a (m+1) x (m+1) system of nonstrictly hy- 
perbolic conservation laws which models multicomponent polymer flooding in enhanced oil-recovery 
process in two dimensions. In the presence of gravity the flux functions need not be monotone and 
hence the exact Riemann problem is complicated and computationally expensive. To overcome this 
difficulty, we use the idea of discontinuous flux to reduce the coupled system into uncoupled system of 
scalar conservation laws with discontinuous coefficients. High order accurate scheme is constructed by 
introducing slope limiter in space variable and a strong stability preserving Runge-Kutta scheme in the 
time variable. The performance of the numerical scheme is presented in various situations by choosing 
a heavily heterogeneous hard rock type medium. Also the significance of dissolving multiple polymers 
in aqueous phase is presented. 



Simulation of two phase flow in porous media plays a key role in many engineering areas such as oil- 
recovery fSHIIESl^ environmental remediation Q and water management in polymer electrolyte fuels 
cells 1 151. We are interested in multi dimensional simulation of two phase flow in heterogeneous porous 
media arising in enhanced oil-recovery. It involves simultaneous flow of two immiscible phases (the 
aqueous phase and the oil phase) in a heterogeneous porous medium. We have assumed that m chemical 
components are dissolved in the aqueous phase. These components could, for example, be different 
polymers that all have different influence on the flow properties. We propose a high order finite volume 
scheme for the numerical simulation of Buckley-Leverett model with multicomponent polymer flooding 
by using the idea of discontinuous numerical flux developed in ||4]|3]- For simplicity we let il = [0, 1] x 
[0, 1] denote the two dimensional reservoir. Let s G [0, 1] denote the saturation of aqueous phase and 
c = (ci, C2, c,n) e [0, Co]™ denote the concentration of the polymers dissolved in the aqueous phase, 
where cq is some non negative real number Then in the absence of capillary pressure the governing 
equations form a (m + 1) x (m + 1) system of hyperbolic conservation laws [24. ,25 J given by 



where {x, t) e fix (0, oo), ai : [0, 1] — > M are given smooth functions and the flux F : [0, 1] x [0, Cq]™ x 

^ M2 js gjygjj by = {Fi,F2), 



1 Introduction 



St +\/ ■ F{s,Ci,C2, C,n,x) = 

{sci + ai{ci))t +V ■ {ciF{s,ci,C2, Cm,x)) = 0, ; = l,2,....,m 



(1) 



Fi{s,c,x) = vi{x)f{s,c 



f{s,c) 



(2) 



A^(s,c) + Ao(s) 
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F2{s,c,x) = [v2{x) - {pw ~ Po)gXo{s,c)K{x)]f{s,c). 



(3) 



Here , po are the densities of water and oil, g is the acceleration due to gravity. The quantities and 
Ao are the mobilities of the water and oil phase respectively and v — (t'l, ^2) G is the total velocity 
given by Darcy law It 161 . 



v = - (^(A„ + \o)K{x)^^, {X^ + Xo)K{x)^ + iX^p^ + XoPo)gK{x)^ (4) 

where K : fl [0, 00) is the permeability of the rock which can be discontinuous in x and p : — > M 
is the pressure. If we assume incompressibility of the flow and if there are no sources, then the velocity 
is governed by 

V-z; = in ^2 (5) 

with some suitable boundary conditions for pressure on dH,. For instance in the inlet part of the boundary, 
water is pumped in at high pressure p ~ pi while a lower pressure p — po is maintained on outlet, see 



Fig 14 On the remaining part of the boundary, the normal velocity is set to zero, which gives a Neumann 
boundary condition on pressure. Equations ([T]l and (|5]l form a system of coupled algebraic -differential 
equations and there is no time derivative involved in equation (|5]l. A commonly used model for the 
mobilities are 

XUs, c) = A„(,s) = ^ >- (6) 

where p^, Po are the viscosities of water and oil and = Pw{c) which is increasing in each of its 
variable The term a; in ([T]i models the adsorption of the component I on the porous medium. 

In the absence of polymer flooding or equivalently if the flux function is independent of c, then this 
problem ([T]) reduces to scalar equation. In lIZTl by using a fast marching method and in ||29l by using semi- 
Godunov scheme method the problem is studied in the absence of polymer Also in [TTj two-phase flow 
problems are studied by using gradient schemes. It is well known that in the heterogeneous media, that 
is when the permeability K{x) is discontinuous , fingering instability fT2l will develop and which results 
in an inefficient oil-recovery. For example see FigfT7|a). As the concentration c increases, viscosity of 



water increases and the fingering effects reduces which leads to an efficient oil-recovery see Fig 17 b). In 
the presence of the concentration c the system ([T]l becomes coupled and non-strictly hyperbolic . When 
the concentration c is smooth, existence and uniqueness theory is established in |37 | but we deal here 
with the case when c need not be smooth. For this system, developing a Godunov type upwind schemes 
are difficult as it needs a solution of Riemann problems. Most often numerical methods requires the 
calculation of eigenvalues and eigenvectors of the Jacobian matrix of the system. Here by using the 
idea of discontinuous flux we reduce the system to an uncoupled scalar equations with discontinuous 
coefficients. Next we study each scalar equation by using the idea of discontinuous flux. This approach 
does not require detailed information about the eigenstructure of the full system. Also in \29\, the idea 
of discontinuous flux is used to study a coupled system arising in three-phase flows in porous media and 
shown its successfulness. Scalar conservation laws with discontinuous flux have been studied by many 
authors |l2l|9l[l0l[ni[l3l[Il[l8l|22l|26l|32l. In particular, in |3| a Godunov type finite volume scheme 
is proposed and convergence to a proper entropy solution is proved, provided the flux functions satisfies 
certain conditions like in ^j2] In one dimensional case for a (2 x 2) system this problem was studied 
in |4| and there proposed a finite volume scheme and named numerical flux as DFLU. This DFLU flux 
works even in cases where the upstream mobility gives an entropy violating solution 1,34]. Here we are 
extending DFLU to a multi dimensional case with high order accuracy. The difficulties of developing an 
upwind type numerical schemes in a highly heterogeneous media in the presence of gravity attracts the 
importance of the proposed work. 
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The paper is organized as follows. From f|2]the idea of discontinuous flux for one dimensional prob- 
lem is briefly explained and also numerical experiments for high order schemes are performed to show 
their efficiency. In §3 two dimensional problem is introduced and the idea of the one dimensional dis- 
continuous flux is extended. Also high order accurate scheme is constructed by introducing slope limiter 
in space variable and a strong stability preserving Runge-Kutta scheme in the time variable |21 1. The 
resulting schemes are shown to respect a maximum principle. Also two dimensional numerical results in 
various situation are shown for a quarter five-spot geometry. 



2 System of equations in one dimension 

The corresponding (m + 1) x (m + 1) system of equations in one-dimension in the presence of gravity is 
given by 



(SQ +a;(Q))t + ^qF(s,ci,C2,....,c„i,x) = 0, Z = l,2,. 



, m 



(7) 



where t>Q and a; G M, (s, ci, C2, Cm) = (s, c) G [0, 1] x [0, cq]™ and 

F(s, c, x) = [w - {pw - Po)gXo{s, c)K{x)]f{s, c). (8) 

In one dimension the solution v of the equation Q reduces to a constant. We assume that the flux function 
satisfies following conditions: 

1. F{0,ci,C2,...,c^,x) =0, F(l,ci,C2,....,c™,x) = w V 1 = 1, 2, ...,m 

2. The function s F{s, ci, C2, Cm, x) is of convex type i.e, has no local maximum in the interior 
of [0,1] X [0,co]™ seeFig[l] 

3. The adsorption term a; = a;(c/) satisfies /i/(c/) = ^(q) > 0, Vq S [0, 1]. 

The case when v = and F does not change sign is studied in Q. Here we assume v need not be zero and 
allow F to change sign, see Fig. [T] In the absence of gravity, Fg = vfs is non-negative or non-positive 
depending on u > or < 0. Hence F is increasing or decreasing in s accordingly. In the presence of 
gravity Fg becomes, 

Fs = ^TT-rWt^ + (P^ - Po)gK{x){s\^ - (1 - s)A„)] 

which vanishes at s = 0, 1. Depending on the values of v, p^, Po,g, K, there can be a root G (0, 1) 
which makes F non-monotone in s, as shown in Fig|T] If such a root exists, it is a root of the following 
cubic equation 

p-w(c) [pu, - Po)gK 

This cubic equation has one real and two complex roots, the real root is given by 

1 



1 + r 



where 

a = -27r + 27r^ - 27z ~ 5Arz - 27r^z, (i = 2916r^ + . 
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Figure 1: A' = A* 



Since 



Fss (-^^ 1 



6341 ~ s*){pw - Po}gK{x) 
(A^ + Ao)2 



then F attains the maximum(minimum) at s = if > po (pw < Po)- Note that the nature of the 
extremum depends only on the densities and is independent of the polymer concentrations c; and the 
permeability K. 

If F{s^ c, a;) — F{s, c) then the system (j7|i can be put in the matrix form as 

Ut + A{U)U^^Q, U =[s ci C2 ... c^Y , 
where A{U) is the (to + 1) x (m + 1) Jacobian matrix 



/dF dF 

ds dci 

^ 



AiU) 







s+hi 





dF 
dc2 



F 



dF \ 



dc, 





\o ... . 

The eigenvalues of this system are given by 

A" = A(s,c) 



dF 
ds 



A' = X\s,c) = 



Fis,c) 



(s,c) 



Z = 1, 2, .., TO. 



s + hi{ci) ' 

We can observe that for any c — (ci, C2, c„i) G [0, 1] x [0, cq]™ and for some / G {1, 2, to} there 
exist at least one point s* = s* (c) G [0, 1] such that (see Fig{T]i. 

A'(s*,c) = A"(s*,c). 

For this couple {s* ,c), A' — A'*, hence eigenvalues may coincide and the problem is non strictly hyper- 
bolic. The Rankine-Hugoniot condition corresponding to (|7]i is given by 

F{s^,c^,x+) ~ F{s^,c^,x-) = a{s^'~s^) 
c^F(s^',c",x+)-cfF{s^,c^,xY = als^cf + ai{cf')^s^cf-ai{cf)) (9) 

V I = 1,2,.., to. 
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For details see [1231,1251. If = c^( i.e. cf — cf V I — 1,2, ...,m) then second equation reduces to 
the first equation of (|9]). This corresponds to the Rankine-Hugoniot condition for single Buckely-Leverett 
equation M. Now we are interested in the case 7^ c^, i.e. cf 7^ cf for some 1,1 < I < m. If we 
combine the two equations (|9]l then we may write 

(cf - cf^)F{s\ c\x-) = a(cf - cf + <T(a,(cf ) - ai{cf)) (10) 

Define the functions hf by 

/^f(Q)= --^'^ V (11) 



/i;(q) if Q 



H ■ 



Now from Q and ( 10 1, finally we get 



G, (12) 



where h — h\ (cf^) (= h\ (cf) V l) . Thus the Rankine-Hugoniot condition reduces to ( 12 1. This gives 
an idea how to obtain a weak solution of the Riemann problem to (|7]l. 

2.1 Riemann problem 

For simplicity we restrict our study to the case when to = 2 in equation ([t]), i.e c = (ci, C2). Also we 
assume that F{s, c, x) = F{s, c). Consider the Riemann problem associated to the system (j7| with the 
initial condition 

^''~[sR if x>0 ' (cf,cf) if x>0. 



Solution to (|7]l and ( 1 3 1 is constructed by connecting states so that it should satisfies the Rankine-Hugoniot 



condition. There are two families of waves that arise in the solution of the Riemann problem referred to 
as s and c waves, s waves consists of rarefaction and shocks (or contact discontinuity) across which s 
changes continuously and discontinuously respectively, but across which both ci and C2 remain constant, 
c waves consists solely of contact discontinuity across which both s and ci , C2 changes such that 
remains constant in the sense of (T2I1. For different choices of cl and cr, the possible shapes of F{s, cl) 
and F{s, cr) are shown in FigK We restrict to the case when cl > CR{i.e., cf > cf', I — 1, 2). When 
Cl > Cr. the fux functions s — >■ F{s,c^) and s — >^ F{s,c^) are one of the shapes given in FigOl To 
explain the Riemann problem, for simplicity we consider the shape of the flux functions as in FigB] 



Case 1: < s* 

Draw a line through the points {—h, 0) and (s* , F{s^ ,0^,02)). This intersects the curve F{s, cf, cf" 
at the point s, where Fs{s, cf , cf ) > 0. We divide this in two subcases. 

Case la: > s 

(a) Connect (s^, cf , cf ) to (s*, cf , cf ) by a s-rarefaction wave (see Fig[3^). 

(b) Connect (s* , cf , cf ) to (s, cf , cf ) by a c-wave with speed (see Figj3^). 

^c=^^^^^#^=F.(.*,cf,cf) 
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Figure 2: Possible shapes of flux functions for different choices of and cjzj. 

(c) Connect (s, cf , c|^) to (s^, cf, cf") by a s-rarefaction wave (see Figjs^). For example if 
F(s, cf , C2 ) and F(s, cf , c|^) are strictly convex functions then the corresponding solution of the 
Riemann problem is given by (see Fig|3]3) 



(s^c^c^) 

((F,)-i(f,cf,C2^),cf,ci') 



if a; < CTgi, 

if (Jst < X < (Jet, 

(s, cf,C^) if (Jet < X < (Tit, 

{{Fs)-\l,cf,4),cf,4) if (Jit<x<a2t, 

„R Ji 



(s«,cf,c«) 



if X > a2t. 



F(s ,c/*,c^ }^ 





'(■y'^cf.cf ) 



(s']cf,C2 ) 



(a) 



(b) 



Figure 3: Solution of the Riemann problem ( 13 1 with s < s* and s > s 



• Case lb: < s 

(a) Connect (s^, cf , C2 ) to (s*, cf , C2 ) by a s-rarefaction wave (see Fig|4^). 

(b) Connect (s* , , cf" ) to (s, cf , c^) by a c-wave with speed (see Fig|4^). 

_ F(s,cf,c«) 



s + ft. 



^s('S J C]^ , C2 ) 



(c) Connect (s, cf, c|*') to (s^, cf , c^) by a s - shock wave with speed (see Fig|4^). 

^ R 

In the case of convex fluxes we can write the solution of the Riemann problem as (see FigHb) 



(s^c^',c|') 



if X < (lit, 



(s(x,t),ci(a;,i),C2(x,t)) 



((F«)-i(f, 0^,4), 0^,0^) if ait<x<act, 



(s, cf^, C2) 
(5«,cf,cf) 



if iJct < X < (Tst, 

if a; > a^t. 




Js.clcf) 



(i.cy.cjj " (s')cf,c^) 



(a) 



(b) 



Figure 4: Solution of the Riemann problem ( 13 1 with < s* and < s 



• Case 2: > s* 

Draw a line joining the points {—h, 0) and (s^, F{s^, cf , C2 )). Let (s, cf , cf^) be the point where 
this line meets the curve F{s, cf , cf), where Fs{s, cf , c|^) > 0. Consider the following subcases. 

• Case 2a: > s 

(a) Connect (s^, cf , C2 ) to (s, cf , c^) by a c- shock wave with speed (see Figjs^). 

_ F(s^,cf,cf)-i^(5,cf,C2«) 



(b) Connect (s, cf , cf-) to (s^, cf , cf^) by a s- rarefaction wave (see Figjs^). 

In the case of convex flux the solution of the Riemann problem is given by (see FigjS]?) 



{s{x,t),Ci{x,t),C2{x,t)) = 



(«^cf,4) 

(s,cf,cf) 



if a; < act, 

if act < X < ait, 



((F,)-i(f,cf,cf),cf,c«) if ait<x<a2t. 



(s«,cf,C2«) 



if X > a2t. 
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(a) 



(b) 



Figure 5: Solution of the Riemann problem ( 13 i with > s* and > s 



Case 2b: < s 

(a) Connect (s^, cf , C2 ) to (s, cf", c|^) by a c- shock wave with speed (see Figj6^). 

" L = 

(b) Connect (s, cf, Cj-) to (s''^, cf , cf) by a s-shock wave with speed (see Fig|6^). 

f(g,ef,c|^)-F(g^,cf,c|^) 
s — s" 

In the case of convex flux the solution of the Riemann problem is given by (see Fig|6]3) 

(5-^,0^,4) if a; < aat, 
{s{x,t),ci{x,t),C2{x,t)) ^ (s, cf,cf) if act<x<ast, 

{s^,c^,c^') if x>ast. 




(a) 



(b) 



Figure 6: Solution of the Riemann problem ( 13 1 with > s* and < s 



Remark: When the flux function F{s, c, x) is smooth in s and c and discontinuous in the x variable then 
the construction of Riemann problem is explained in the appendix of [4]. Here also we can construct the 
solution of Riemann problem in a similar way. 
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2.2 Finite volume scheme 



We define the space grid points as x^^i = ih, h > and i € Z and for At > define the time 

discretization points t„ = nAt for all non-negative integer n, and A = The Finite volume scheme 
for the system (j?]) is given by 



n+l n+1 I / n+l\ 



ci^r + ai(cir)-A(G, 



(14) 



C2 



n+1 n+1 



+ «2(c2r') = c^ls^ + a,{c,f) XiG^^ -G2IO 



where the numerical flux F" 1 , G" 1 and G" 1 are associated with the flux functions F(s, c, x) and 
G/(s, c, x) = ciF{s, c, x), I = 1,2 and are functions of the left and right values of the saturation s and 



the concentration c at x 



jpn TP/' n n n n n n ^ \ n /-i f n n n n n n- ^ \ 

Fi+i = ,ci, ,C2, ,s,;+i,ci,+i,C2,+i,a;,+ i), Gl^_^_l = Gi{Si , ci, , C2, , s,;+i, ci,+i, C2,+i, a:;,;+i ). 
The choice of the numerical flux functions F and Gi{l = 1,2) determines the numerical scheme. Once 



we compute s"^^ from the first equation of (14 1 then we recover ci"^^ and 02"^^ from second and third 



equation respectively using an iterative method, like Newton-Raphson method. 
Now we briefly explain the DFLU flux of [4J and Godunov flux. 

2.3 The DFLU numerical flux 

The DFLU flux is an extension of the Godunov scheme that was proposed and analyzed in IS) for scalar 
conservations laws with a flux function discontinuous in space. We define 



Gn 



n pn 

n 



if F" 1 > 

1+5 

if F" 1 < Z = 1,2. 



(15) 



Now the choice of the numerical scheme depends on the choice of F"_^^^r^. To do so we treat c{x, t) in 
F{s, c, x) as a known function which may be discontinuous at the space discretization points and F is 
allowed to be discontinuous in the x variable at the same space discretization points. Therefore on each 
rectangle {xi_i, ^i+i) x {tn, tn+i), we consider the conservation law: 

St+ F{s,ci^,C2l,x)^ = 
with initial condition s{x,0) = s" for x^_i < x < x^_^_i (see Figj7|. The above problem can be 

t — tn+l 





St+Fis^Ci'l, 02^1, x,)^ ^ 


St+F(s,Ci"+i, 02,"+!, 2:^+1)2: = 













Figure 7: The flux functions F{., ci,C2,x) is discontinuous in ci, C2 and x at the discretization points. 
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considered as a conservation law with flux function discontinuous in x for which DFLU flux can be used. 
Then the DFLU flux is given as 

TjDFLU ( „n „ n n n n n \ 

= max{F(max{s7, 6*,"}, ci^, cz?, Xi), -F(min{s^+i, 6l,'Vi}, ci^!^_i, ca^+i, x.+i)}, 
where 9^ argmini^(., cif , C2", x^). 

2.4 The Godunov flux 

The Godunov flux at the grid point x^j^ 1 is calculated by using the solution of the Riemann problem: 

st + F{s,c,x)x = 

{sci+ai{ci))t + {ciF(s,c,x))x = 

{sc2+a2{c2))t + {c2F{s,c,x))^ = (16) 
in the domain 1 , x^^i) x (t„, t„+i), with the initial condition 



(six, t„), ci(x, i„), C2(x, t„)) 



(sf,cif,C2r) if a;<a;i+i 

(sr+i:Ci?+i,C2r+i) if a;>Xi+i. 



The numerical fluxes are given by 

■Fl+i = F{s{x^+l,t),c{x^_^^l,t),x^_^_l), t.n<t< i„+i 

and 

Remark: In general Godunov and DFLU flux may differ, for details see H . 

2.5 The Upstream Mobililty flux 

This flux is designed by petroleum engineers from physical consideration. It is an ad-hoc flux for two- 
phase flow in porous media which corresponds to the approximate solution to the Riemann problem [81. 
To define the upstream mobility flux, assume that the absolute permeability K{x) > and we redefine 
the flux function in ([8]) as 

KX 

F{s,c,x) = — /" [v - {pn, - Po)gKXo{s,c)] (17) 

Now if we take A^, = KX^^ and Ao = KX^,, the flux function becomes 

F{s, c, x) = [v - {pyj - po)gXo{s, c)] 

Now the numerical fluxes are given by 

i (^r, c?, , c^, , , c^'^+i , cj,^ J = [v ip^ p,)gX:] , 



Xt 



( Xe{s'^,c'^^,c^^,Ki) if V - {pe " Pi)gXe > 0, i ^ w,o,i ^ e, 

I A£(sf+;^,cJ-^_^^,c5,.^pXj+i) if V - {pi- p^)gXi <0, i = w,o,i ^ £ 
and 1 (Z = 1, 2) are given as in (jTsj). 



Remark: The Upstream mobility flux works only for the flux function which is of the form as in ( 17 1 
where as DFLU flux can be applied for any flux function which satisfies the assumptions of fj2] 
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2.6 High-order schemes 

In order to develop the second order scheme, we follow the method of lines approach in which space and 
time discretization are performed separately. In the first step, spatial discretization using piecewise linear 
reconstruction is made which leads to a system of ODE which can be written as 



dU 

~dt 



R{U) = 0, U 



sci + ai(ci) 

SC2 + 02(02) 



1 

Aa; 



' 2 
Gij_ 1 



(18) 



The high order accurate fluxes are given by 



G. 



(19) 



The quantities with superscripts L and R denote the reconstructed values of the variables to the left and 
right of the corresponding cell face. For any quantity u, we can define the reconstruction as follows: 



1, 

= + 7:0, 



1 



(20) 



where 



5i = minmod 



u,_i),6I(m,+i - u,) ,6* e [1,2] 



(21) 



Finally the time integration of the ODE (18 1 must be high order accurate in order for the scheme to be high 
order accurate. A third order accurate, strong stability preserving Runge-Kutta scheme due to Shu-Osher 
is given by 



jjn+l 



jjn 

- Aii?(F(o)) 
^C/" + i[l/(i)-Ati?(T/(i))] 

o o 



If the explicit scheme (14 1 is stable in the norm || . ||, i.e., if 

At < Atr = 



U - AtR{U)\\ < \\U\\ , where AU is the CFL restricted time step. 



(22) 



then the above Runge-Kutta scheme is also stable in the same norm under the same time-step restriction 

(cf.cniEii). 



2.7 Maximum principle on saturation 



Let us write 



nR 



nL nR 



I) 



nR 



nR 



„ nR nL nH \ 
C2,_i,C2i+l,C2,^l j 



nR 
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The updated value of the saturation ( 14 1 can be written as 

Where H is Lipschitz continuous in saturation and concentration. Since the slope limiter preserves the 
average value of the solution in each cell, we can express this as 

s^'+i = ^ '"^ - XiFll, - i^ii). (23) 
If we differentiate H with respect to its variables s" we can observe that gfrri? > provided 



A||r^4,l<^. (24) 



Let 

dF2 Fi F2 . 
M = sup{ — , — , —-J-}, 
s OS OS s + hi s + hi 

then the condition ( |24) l reduces to, 

AM < -. (25) 
- 2 

This shows that H is monotone in each of its variable. Using these facts we have the following lemmas. 

Lemma 2.1 Let sq e [0, 1] be the initial data and let {s"} be the corresponding solution calculated by 
the finite volume scheme ( |i4p using DFLU flux along with slope limiter If the CFL given in ( |25p holds 
then 

< s" < 1 \/,iand n. (26) 

Proof: From the property of slope limiter we can observe that whenever < s" < 1 then the recon- 
structed values satisfies 

< s"fi,s"fi < 1 V iandn 
Using this property and the monotonicity of the H , we get 

= i7(0, c") < iJ(s", c") = s^+i < c") = 1 

This proves that 

< s"+^ < 1 V 



2.8 Maximum principle and TVD for concentration 



Theorem 2.2 Let {ci"} , {c2"} be the solution calculated by the finite volume scheme ( 14 ) using DFLU 
flux with slope limiter Under the CFL condition AM < concentration c = (ci, C2) satisfies 

(a) min{cjf_i,Q^Qf+i} < c4+' < max{Qf_i, , Qf+ J V n G Z+, z G Z 1 = 1,2. 

(b) El'^'^'-^'-i'l^El^'"-^'"-!! VneZ+, 1 = 1,2. 
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Proof: From the finite volume scheme ( [T4] t, 

„n+l _ n \( TPn ZPn 



We can express the numerical flux Gij^^ i , G2i^i as (here we suppress the index n for fluxes ) 



where 



f^+i =max{F,+ i,0}, F^, = min{f^,+ 1 , 0} 



We write the scheme (14^ as 



By adding and subtracting the term s"+^q", we get 

where a/(c/^+i)-a,(Qr) =a;(Cr^*)(c;r^'-Cir), for some Cr^Hetweenc4'+i and c;f. By replacing 
s^' - by -A(F,+ i - ) and splitting 1 by (^;+ ^ + J^'i ) we have 

t-r 2 '■I 2 2 2 

+ A(qI^. K^. + c^^k;, - ic^\Fl^ + c4^^ J) = 0. 
By rearranging the terms in the above equation we get 

Note that 

^nL _ ^ nR _ ^ n 



and 5i is the slope limiter given by 



5, = minmod ( e{c{: - c^ti), ^{ci!l+, - c,ti),0{ci:+, - ci^) 



1 

After substituting the values for ci^j^i and ci"]^i the above equation becomes 

(-r^ + aKcr - + a^^+ , | + af- , (i - ^^^^^^^i ^^^ )ioi7+i - o 

+ XF+ 1 (1 - ^ ^)(c;" - q" 1) + XFr 1 - = 0. 
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Now we can write 



ci^' = ci7 - A ^ 5M7 - citi) (27) 

F~ 1 § 



n / I |3N/n n\|/2 ,4\/'n n\ 

= c;? - (qr - citi) + I^r+|(c4Vi - c;?), (29) 

where 
and 

p+ p— 



From the properly of the limiter it is easy to see that 

< — — r < 1 (30) 

-2(Qr+i-Q?)- 

which in turn implies 

Cf_i,D,"+i>0 (31) 
Now we prove the maximum principle for ci{l = 1,2) by considering the following cases. 
Casel: Suppose that hes between ci2-\ and Q -^.^ then 

c{i = eciti + (1 - ^)c,r+i for some 9 e [0, 1] 

and 

Ci7-ciU = {^-0){c,7+i-ci7-i) 
cif+i - c{} = 0{ci^+, - ciU)- 
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Now from ( p9] l we write 

Q^+i = (1 - 0)(qIVi - - - ^)(c;^Vi - 

+ ^,'Vi%;r+i-c;r-i) 

where 



(32) 



Ai =0(l-Z?f^i) + Cf_i(l-0) 

A2-(i-e)(i-qLi) + e^r+- 

t 2 ^ 2 

Note that Ai + A2 = 1, under the CFL condition AM < 5 we have C^_i,D'^_^i < 1 which gives 



A; , A2 > 0. Hence from ( 32 1 the maximum principle (a) follows. 

Casel: Suppose q" does not lies between and q"^]^, then we have 5i — 0. i.e., 



The equation (29 1 can be rewritten as 

„ n+l /I ^ri j-,n A„ " i /ti „ n i n" ^ " 

Q, - (1 - c,^_i - L),+ i)cii + c._iQ,_i + 

Note that 

C" 1 + LI" 1 < 1 under the CFL condition AM < ^ . 

1-5 »+2 - - 2 

This proves the maximum principle(a). 
To prove the TVD property, consider 



'Ji+1 



(sr'+«i(C^))2(Q^i-Q?) 



— (1 



A 



< AM 



2(qIVi - O 2(c4Vi - 



-1 



2(qLi - Q?) 2(q 



= 2AM < 1, 



(33) 



(34) 



under the CFL condition AM < i. From (31 1 and (34i the TVD property(b) follows from the Harten's 
lemma. I 
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Figure 8: Saturation s for first order scheme (left), high order scheme (right) at time < = 1, mesh size 



" 100 • 



c 

of s. 



Remark: Note that saturation s need not be of total variation bounded because of / = /(s, c) and 
c{x, t) is discontinuous (see fTl). The singular mapping technique as in |3| to prove the convergence 



J looks very difficult to apply. However by using the method of compensated compactness, Kaksen, 
Mishra, Risebro I 



showed the convergence of approximated solution in the case of a triangular system. 
By using their results in the case of a single component polymer [m = 1) under suitable assumptions, in 
Pi convergence analysis of the saturation is studied. 



2.9 Numerical results 

Here we have chosen the flux function for the above system of equations (j7]l with v ^ Q.2 , K = 1 , 
= o.s+e^+e^ , Ao = (1 - s)^, p^g = 2 and pog = 1. The adsorption term is given by a;(c;) = 
1 + 0.5c/ {l = 1, 2). In the numerical experiment the initial data is chosen so that the flux function F is 
allowed to change the sign, equivalently eigenvalues X\l — 1,2) of the system (j?]) allowed to change the 
sign. For this purpose the initial data is chosen as 



(s(a;,0),ci(a;,0),C2(x, 0)) = 



(0.1,1,0.6) if x<OA 
(1.0,0,0) if x>OA. 



Numerical experiments are done for DFLU flux. Upstream mobility flux and compared with Godunov 
flux. In these experiments data are chosen so that DFLU flux differ from Godunov flux. The performance 
of the DFLU flux is as good as the Godunov flux. High order accurate schemes corresponding to DFLU 
and Godunov are constructed by introducing slope limiter in space variable and a strong stability preserv- 



ing Runge-Kutta scheme in the time variable, a comparison with first order scheme is shown in Fig 8]9 
For first order and high order scheme it clearly shows that the DFLU flux is as good as the Godunov 
flux. Note that Godunov flux requires the solution of the Riemann problem of a system where as DFLU 
flux requires the solution of the Riemann problem of a scalar equation. Fig |10|l l|12| shows that numerical 
solution computed by DFLU is as good as Godunov and converges faster than Upstream mobility scheme. 

For high order scheme the error , order of accuracy a for DFLU, Godunov and Upstream mobility 
schemes are given in the table 1 . The order of accuracy a is calculated as foUows: 

lrL(ci /c2 ) 

ei = ||s - s/iJlii with /ii = /i, e2 = lis - s/i,,|lii with /i2 = /i/2, a= — — - — . 

In 2 
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Figure 9: Concentration ci for first order scheme (left), high order scheme (right) at time t — 1 (right), 
mesh size h — . 
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Figure 10: Saturation s at time t ~ \ with mesh size h — (left) and h — (right), with high order 



accuracy. 
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Figure 1 1 : Concentration ci at time t = \ with mesh size h ~ (left) and h = ^ (right), with high 
order accuracy. 
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h 


DFLU 


a 


GODUNOV 


a 


UPSTREAM 


a 


1/50 


4.2336x10-^ 




4.8839x10-^ 




6.3189x10-2 




1/100 


2.4366x10-^ 


0.7970 


2.7735x10-^ 


0.8163 


3.6055x10-2 


0.8095 


1/200 


1.3605x10-2 


0.8407 


1.5268x10-2 


0.8612 


1.9805x10-2 


0.8643 


1/400 


6.2334x10"^ 


1.1260 


6.9589x10-3 


1.133 


9.2108x10-3 


1.1045 


1/800 


2.2233x10-3 


1.4873 


2.4398x10-3 


1.5121 


3.3674x10-3 


1.4517 





DFLU 




GODUNOV 




UPSTREAM 




h 




a 


l|ci-c?|Ui 


a 


I|C1 - c'IWli 


a 


1/50 


3.3257x10-2 




3.9971x10-2 




5.0529x10-2 




1/100 


2.2303x10-2 


0.5764 


2.5938x10-2 


0.6239 


3.4946x10-2 


0.5319 


1/200 


1.2304x10-2 


0.8582 


1.4014x10-2 


0.8881 


1.874x10-2 


0.899 


1/400 


4.8878x10-3 


1.3318 


5.4714x10-3 


1.3569 


7.9071x10-3 


1.2449 


1/800 


1.6586x10-3 


1.5592 


1.8413x10-3 


1.5712 


2.8197x10-3 


1.4876 






DFLU 




GODUNOV 




UPSTREAM 




h 


I|C2 -4\\l^ 


a 


I|C2 - C2IILI 


a 


\\c2-4\\l^ 


a 


1/50 


1.9954x10-2 




2.3983x10-2 




3.0318x10-2 




1/100 


1.3382x10-2 


0.5764 


1.5563x10-2 


0.6239 


2.0968x10-2 


0.5319 


1/200 


7.3821x10-3 


0.8581 


8.4086x10-3 


0.8881 


1.1244x10-2 


0.899 


1/400 


2.9327x10-3 


1.3318 


3.2829x10-3 


1.3569 


4.7455x10-3 


1.2445 


1/800 


9.9518x10-"* 


1.5592 


1.10479x10-3 


1.5712 


1.6924x10-3 


1.4875 



Table 1: error for saturation s and concentrations ci and C2- 
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Figure 12: Concentration C2 at time t ~ 1 with mesh size h = ^ (left) and h — (right), with high 
order accuracy. 
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Note that as /i ^ 0, a in DFLU is better than a in Upstream and more close to a in GODUNOV. Here 
the exact solutions s,c — (ci, C2) are computed from Godunov scheme for very small values of h and At 
with j^M = 0.5. 

3 2-D model 

In this section we are extending the numerical schemes explained in § 2 for one dimension to a multi 
dimensional space. For simplicity we explain only in two dimensions and higher dimension can be handled 
in a similar way. In dimension two the equation ([T]) can be rewritten as 

., + f|i(.,c,x) + f|i(.,c,x) = 
{sc,+a,{c,))t + ^is,c,x) + ^{s,c,x) = (35) 

where (x, t) e ^ x (0, 00), x = {xi,X2) and the flux i^i, F2 : [0, 1] x [0, co]^ x 17 — >• M are given by 

F,{s,c,x) = v,{x)f{s,c), f{s,c)^ /"("'^l (36) 

Xyj{s,c) + Ao(s) 

F2(s, c, x) = [v2{x) - {pw - Po)gK{s, c)K{x)]f{s, c) (37) 

To compute Fi and F2 we need the velocity component v — (f 1, U2). This velocity (pressure) is governed 
by the incompressibility of the flow: 

V • w = in n (38) 

with some suitable boundary condition for velocity (pressure) on dil as explained in § 1 . 
Basic numerical approach for finite volume method is outlined in the following algorithm: 

1. Set time step n — and initialize s°. — (c°, C2). 

2. Assume s" and c" ~ (c" , C2 ) are known att — tn- 

3. Solve for the pressure from (|4|i and (jSj. 

4. Compute velocity from (Hh. 



5. Chose time step At" so that CFL condition is satisfied see { 3.5 . 

6. Update saturation and concentration at t = tn+i level by 



s 



ci ^ + ai{c^ ) = s ci +ai(ci) — At V • (ci _f (s ,c ,w jj 



s"+ic2"+i+ai(c^+i) = s"c2"+a2(4'+')-At"V-(c2"f^(s",c",t;")) 
7. Set n = n + 1 and Go to step 2. 
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(i4.H) 



Figure 13: Definition of cell Qi j by grid points 



3.1 Discretization of the domain Q = [0, 1] x [0, 1] 

Consider the Cartesian grid obtained by taking the cross product of the one-dimensional partitions {xi, i = 
1, . . . ,nx} and {yj, j = 1, . . . ,ny} with xi — yi ~ and Xn^ = yuy = 1- We also introduce 
one layer of grid points on all four sides of fi which will be referred to as ghost points. Thus the 
grid point indices range over < i < n^, + 1 and < j < Uy + 1. The grid defines the cell 

QiJ — [^j-ij^^i+i] X 1 , y,;^i], see Fig 13 for < i < rix and < j < Uy. The number of 



true cells where the solution is supposed to be computed in the domain VL is iic = {rix — 1) x {riy — 1) 
(excluding the ghost cells). 



3.2 Numerical approximation for the pressure 

Define /i := (A^ + Xo)K and 9 := {XwRw + XoPo)gK. Integrating equation (jsj over cell j and using 
the divergence theorem, we obtain the finite volume approximation 



where the velocity at the cell face is given by 



dp 

v.^U = -^Yx 



dp 
dy 







We approximate these as follows: 
• Along the x-direction 



Pi. 



^j.- da; ; 



dp 

dx 



dx /i 
1 



dp 

'' dx 



1 



1 



-da; 



Ax 



This leads to the following approximation for the velocity flux 



Pi+i,j - Pij 
Ax 



1 



1 



(39) 



(40) 



• Along the y-direction 

Pz.j + l -Pi.j 



'V, 



dp 
dy 
Ay ( 1 



dy ; 



1 



Ay 
2 



Vj 
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Hence we get the approximation 



i - + i ^ ^.j+i (41) 



where 



(42) 



The velocity on the inlet boundary is computed as 



Pi,j ~Pi 1 1/^1,1 



with similar expressions for the other inlet/outlet parts of the boundary. On the rest of the boundary, the 
normal velocity is zero which is equivalent to saying that flux is zero. The system of equations (|39]) for 
the pressure can be put in the form 

Ap^h (44) 
where A E ]R"<:>^"<: and b G M."". This matrix equation is solved using conjugate the gradient method. 

3.3 Finite volume scheme 

By integrating equations in ([T]) over the cell Qij, we obtain the following finite volume approximations 

s"+i=s'\---^^^|[F"i .-i^" 1 .lAu+fF".^! -F". i]Ax\, (45) 



AxAy 

+ lic^F)l^^^-{c,F)l^_^]Ax}, (46) 



+ [{c,F)l^^,^~{c,F)l^_^]Ax]. (47) 

Here we introduce the DFLU numerical flux for two dimensional finite volume scheme by using the idea 
explained in §2. The corresponding numerical fluxes are given by 

where c"^- = (ci"^ , 02"^ ) and 

{e\),^j = argmin Fi{.,clj,Kij) and (6*?.^^^ = ai-gmin F2(., c.^^^, 



Q, 1 . ifF" 1 . > 

cu+1 iF"]! ■ ifF" 1 . < ? = 1,2. 
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3.4 High-order scheme 

In order to develop the second order scheme, we follow the method of lines approach in which space and 
time discretization are performed separately. In the first step, spatial discretization using piecewise linear 
reconstruction is made which leads to a system of ODE which can be written as 



where 



dU 

'dt 



1 

AxAy 



F{s, 



1 ^,;4- 



+ RiU) -- 


= 0, 


U = 


SCl 
_SC2 


S 

■f ai(ci) 
-f a2(c2)_ 




[Gii+i J 


- Gij 

-G2, 


-hJ^y + 






given by 
















1 ,-,ci 

2 'J 


R 







(48) 



(49) 



Gi 



if P 



+ k.3 



if P 



> 

<o ; = i,2, 



(50) 



and similar expression for ^j^i. The quantities with superscripts L and R denote the reconstructed 
values of the variables to the left and right of the cell face. For any quantity u, we can define the recon- 
struction in a;— direction as follows: 



where 



5=" . 



minmod 



-(5^ 
2 



(51) 



, ^e[i,2]. (52) 



Similarly in the y— direction we can define 1 and 

3.5 Stability results 



Let us write 



n,L n,i? n,L n,H n,L n,R n,L n,R \ 
S . 1 .,5.1 ..S. • 1,5.. ii) 

~ - '^jj' «J-5 «J-3 + 5 + 



*- 2 J 



/-n\16 / Ti,^ ^T'.rt n.iv n,/t, 

(C, = ^.,Ci^_, ^.,Q^_^, ^.,Ci^^i 



-,- n,R. N 



The updated value of the saturation (45 1 can be written as 

Where is Lipschitz continuous in saturation and concentration with the property 
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Since the slope limiter preserves the average value of the solution in each cell, we can express this as 

n,L , n,R n.L , n,_R 

- 4 + 4 Ax^^^+hi " ^-i.^-^ " Ay^^^^+I " ^^^-i^' ^^^^ 

If we differentiate H with respect to its variables s" we can observe that > provided 

Let 



Af = sup{— — , — — , , , , , , I, 
s as as s + hi s + hi 



then the condition (|54| reduces to, 



nmx{A="M, A^Af} < ^ where = A^ = (55) 
4 Ax Ay 

This shows that H is monotone in each of its variable. Using these facts we have the following lemmas. 

Lemma 3.1 Let sq G [0, 1] be the initial data and let {sf j } be the corresponding solution calculated by 
the finite volume scheme \45\ using DFLU flux along with slope limiter If the CFL given in \55\ holds 
then 

< < 1 V, i, j and n. (56) 

Proof: From the property of slope limiter we can observe that whenever < sfj < 1 then the recon- 
structed values satisfies 

< s";^i'^, s"'.^'f < 1 V i, 7 and n. 
Using this property and the monotonicity of the H , we get 

= H{0, c", Wi± 1 J, w.ji 1 , K,±ij,K,^j,K,^j±i) 

< iJ(s",c",W,±l^^-,W,^j±l,K,±ij,if,j,if,,,±l) = s^+1 

This proves that 



Now we prove the lemma that gives the maximum principle for the concentration. 



Lemma 3.2 Let {c"^ j},{c2 ^ j} be the solution calculated by the finite volume scheme (46 \ and (47\ by 
using DFLU fiux with slope limiter. Under the CFL condition ( 55 i concentration c — {ci , c^) satisfies the 
followig maximum principle 

(a) min{cr,,,.,cr,±i,,.,cr,,,.±i} < ci-+^ < max^,,,., cr,±i,^., cr,,,.±J, 

y neZ+, ieZ,l = l,2. 
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Proof: We can express the high order numerical fluxes Gl^_^_l j, Gi^ j+ 1 = 1, 2) in the finite volume 
scheme (|46]) and as 



1 , =cr^, ^+1 ■+c?'\ F. 



where 



K^,^^ =ma^{^;+.,^,0}, K;, ^ =min{^;+.,^,0} 



We write the scheme (46 1 and (47 1 = 1, 2) as 



+ l„n+l , /^ri+lN „n „n „ / „n \ , \2:/_-L p+ i „B. jp— „L tp+ „R 

3 ~«/K,) + A (q^^i 



By adding and subtracting the terms s"^^c"^ we get 

wherea;(Q"+^)-a;(c/^j) = a;(Crj^" )(c/"+^-Qfj), for some Cf/' betweenQ"+^ and c/^^. Replacing 
^"^ - by -A-(^;+i, - . ,) - Anj;.,+| - . ) and splitting . by , _^ + , 
(similarly for F^ ) and by rearranging the terms we have 

+ A"F.+ 1 .(c^ . - , ) + A^Fr 1 .(cr . - cr"^, ) 
+ A''^;:,.. (ci:,, - cr;_ J + A^^-,_. - cr;;_ j = o. 

Note that 

nL _ n I nR _ n 1,3 „nL _ n , i,3 „nR _ n i,3 

H.^l . — Hi , + r, ^ H. 1 . — hi i o ' ~ o ' . 1 — hi i o 

1+^.3 2 i-^.j '-J 2 '.3 + 5 "'^ 2 2 

and 5f J and ^ are the slope limiter given by 

Sl^ = minmod (^9{cl^ - cl_J, \{cl^^^^ - cl_J,e{c^^^^^^ - cl j^ , 

61^ = minmod (^0(ci: ^ - c^.^.^ J, ^(c^^,^,, - cr.,^_ J, 0(cr^.^.^, - . 
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After substituting the values of , c, , c, , c, the above equation becomes 

'«±lj 'i±l,j 'ij±l 'i,3±l 

=cl . - A" ^(cf -cl , .) 

(s^j + O'KCj ))(cr,,^- - cp^_,, .) 

Kr+«Kc::;^))(ci:,,.-ci:„) ' ''^ 

- A^' ^(cr -cf _ ) 

1 _ 
1 5^ 

(C+'+«KC'))K,,,-cr,,) 2 



Now we write it as 



2 

j + i 

-c". i(cr -cv )+£'"._^i(cr -cr ), 



where 



C" 1 . = 1 . + 1 ., D" 1 . = 1 . . 

«-5J '+5 J '+5 J 

*J-5 «J-5 «J-2 «>J + 5 «J+5 «J + 5 
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and 



Kr+«KC^))K,-cr,_,,) 2 ' 

1 ■ 

1 -^A- ''''' 1 (1- . 



■ '^~l~2\\/' 77, 71 \ ' 

F~- , 1 

^3 _ \y "'J~2 /-r 'I 

1 <5^. 

(^"r+«;(c;:;b)(ci:,^,,-ci:,p ^ • 



From the property of the limiter it is easy to see that 

- 2(c? - ) ' 2(cr - cr ) ~ ' 
which in turn impUes that 

- 

Now we prove the maximum principle through following cases. 
Case 1: Suppose that 

(a) cI* lies between . and cf and 

(b) cV lies between c? , and , then 

cj*^ ^. = 6l"'cJ^_^ , + (1 - ^'')cr,+i , for some ^ e [0, 1] and 
cr.,, = + (1 - 0'')cr,, .^^ for some 0^ G [0, 1]. 

Now 

cr,,-c, = (i-n(ci:^,,,-cL.,,), 
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By writing c" , as ^ (c" . + c" . ) and substituting the values from ( 58 1 and ( 59 1 the equation ( 57 i becomes 



K,,, - cL,,) + i(^^cr^,,_, + (1 - ncr.,+j 



where 



(1 



Ai 

A2-(i-n(^-q^,,) 



D 



Note that Ai + A2 + A3 + A4 = 1 and with the CFL condition (|55]l we have 



C" 1 .,C". 1,1?" 1 ,£'"^1 < I 

»-5J' »J-5' *J+5 - 2 



which gives Ai, A2, A3, A4 > 0. Hence from (6O1 the maximum principle follows. 
Case2: Suppose that 

(a) c? does not lie between c? , and c? and 

(b) cV" does not lie between c? , and c? , then we have Sf 



0. i.e.. 



a. , 1 . , O ■ ■ ] 



. 1 and Z?" 1 



The equation (57 1 can be written as 



'i.j ^ »-3J «+3J «'J-5 'J+s' 



Note that 



C" 1 ■+D'"i - 



under the CFL condition 



(60) 



This proves the maximum principle. Other cases can be handled in a similar way and the maximum 
principle can be shown. 



3.6 Numerical experiments 

For numerical simulation we have chosen an example of the quarter five-spot problem in the domain 
[0, 1] X [0, 1] . To show the effect of gravity numerical experiments are performed in the presence of gravity 
as well as in the absence of gravity. Also to study the polymer flooding effect numerical experiments are 
performed for various concentration of the polymers. The behavior of water saturation is studied when 
the polymers are injected with different concentrations. The flux function F = (^1,^2) takes the same 



form as in equation ( 36 1 and (BTJl with 



A,, 



-, Ao = (1 - sf, pnig = 2 and pog = 1, ai{ci) = 1 + 0.5c/ (/ = 1,2) 



and velocity v across the grid point is calculated by using (|39 
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P=PI 

(0,0) 



(1.0) 



Figure 14: 
conditions 



Reservoir domain and boundary 



Figure 15: 
boundary 



Pumping water through the inlet 



3.7 Initial and boundary conditions 

The simulations are performed in a computational domain O = [0, 1] x [0, 1] for t S [0, 1]. The initial 
condition is s{x, 0) = 0, i.e. In the inlet part of the boundary we pump water with a pressure p ^ pi and 
we keep the outlet part of the boundary with a pressure p = po{pi > Po)^ on the remaining part of the 



boundary normal velocity is set to zero. The initial inlet saturation is shown in Fig 14 and 15 



3.8 Permeability of the porous media 

We consider a heterogeneous porous medium with an absolute permeability K{x). In order to illustrate 
the robustness of the proposed numerical scheme we consider the two model porous media. The first test 
case corresponds to a heterogeneous medium with a continuous random permeability given by 

N 

K{x) = min{max{V<?,(x),0.5},1.5} (61) 



and 



i=0 



where Xi are TV randomly chosen locations inside the domain. Here we have taken N ~ 100. The second 
test case corresponds to a heavily heterogeneous medium with hard rocks and the permeability is given 
by choosing N random locations Xi and 

K{x)-l^'^^ ifx e B(j;j,0.0015) forsomei e {1,2, ,,7V} ^^^^ 
I 1 eslewhere 



The permeability fields for these two test cases are shown in Fig 16 



Experiment 1: Simulations in this experiment was performed using the spatial permeability distri- 



bution given in (61 1, shown in Fig 16 a). The viscosity of water is given by yUu,(ci, C2) = 0.5 + ci + C2. 
We inject water through the inlet boundary with an inlet pressure pj = 8 and inlet concentration ci — 
and C2 = 0. This is the case of without polymer As expected it produces fingering effects, consequently 
when the water front touches the outlet boundary, a large amount of oil is stuck in the remaining portion 
of the domain, which reduces the efficiency of oil-recovery, this is shown in FigfT7|a). To avoid this 
instability polymer is dissolved with water and injected through the inlet wall. In the presence of polymer 
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(a) 



(b) 



Figure 16: (a) Permeability fields for (61 1, (b) Permeability fields for (62l . 
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(a) 



(b) 



Figure 18: (a) Saturation s in experiment 2, in the abscence of polymer (b) Saturation s in experiment 2, 
in the presence of polymer. 



say ci — 7 and C2 = the fingering instability almost disappears and the amount of oil produced at the 



recovery well (outlet boundary) is increased. This is shown in Fig 17 b) 



Experiment 2: The permeability fields is chosen as in expression ( 62 1 with the presence of gravity, 



shown in FigfT6|b). Which corresponds to a heavily heterogeneous media with hard rocks. Here we have 
taken viscosity of water as /i^(ci, C2) — 0.5 + ci + C2- The result obtained in Fig 18 a) corresponds to 



the saturation profile with the inlet concentration ci = and C2 = 0. The result obtained in FigfT8|b) 
corresponds to the saturation profile with inlet concentration ci = 5 and C2 = 3. A consistent behavior 
of the saturation profile shows that our proposed scheme works well with varying spatial discontinuity in 
the media. 

Experiment 3: This experiment is mainly to study the effect of gravity in saturation profile. This 
experiment is performed using spatial permeability distributions given in ( [62] ), shown in FigfT6|b). Vis- 
cosity takes the form /i^(ci, C2) = 0.5 + ci + C2. We chose the inlet concentrations as ci = 7 and C2 ~ 0. 



The expression involving gravity term is considered along the y direction (see eqn (37 1). The resulting 



figures are shown in Fig 19 a) with the absence of gravity and FigfT9|b) with presence of gravity. Observe 
that presence of gravity significantly effects the saturation profile. 

Experiment 4: This experiment is to study the effect of adding more than one polymer with different 
concentrations. In this model we have taken fiw{ci,C2) — 0.5 + ^/ci + y/c2 and permeability field is 
chosen as in ( [62] i. Figure 20 a) corresponds to the case with concentrations ci 49, C2 = 0. In figure 
20 b) we have taken the concentrations to be ci = 25, C2 = 24. Observe that the total amount (ci + C2) of 



injected concentrations in both the case are the same. But in the second case by adding two concentrations 
the sweeping profile of water saturation is improved considerably. This is reflected in fig[20{b). It is clear 
from this fact that by adding multiple polymers and by taking a suitable viscosity /iu,(ci, C2) it may be 
possible to maximize the oil-recovery. 
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(a) 



(b) 



Figure 19: (a) With out the effect of gravity (b) With the effect of gravity 




(a) (b) 



Figure 20: (a) Single component (b) Multicomponent 

4 Conclusion. 

A high resolution finite volume scheme is developed to study the two-phase flow in porous media by using 
the idea of discontinuous flux. The idea of discontinuous flux helps to reduce the system to an uncoupled 
scalar equation with discontinuous coefficients. Discontinuous flux uses the solution of the Riemann 
problem of the scalar equation where as the Godunov flux needs solution of the Riemann problem of 
the coupled system which is difficult to construct especially in the presence of gravity, heterogeneity and 
multiple components. The results obtained from the idea of discontinuous flux agrees well with the results 
obtained from the Godunov flux. The two-phase flow is studied in the presence as well as in the absence 
of gravity. It is shown that the presence of gravity affects the saturation profile. Also the efficiency of 
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the numerical method is demonstrated by performing numerical simulations corresponding to two-phase 
flow in heterogeneous media. 
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